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Abstract 

This paper presents a hierarchical Bayesian approach to the estimation of com- 
ponents' rehability (survival) using a Weibull model for each of them. The proposed 
method can be used to estimation with general survival censored data, because the 
estimation of a component's reliability in a series (parallel) system is equivalent to 
I the estimation of its survival function with right- (left-) censored data. Besides the 

■ Weibull parametric model for reliability data, independent gamma distributions are 

. considered at the first hierarchical level for the Weibull parameters and independent 

uniform distributions over the real line as priors for the parameters of the gammas. 
In order to evaluate the model, an example and a simulation study are discussed. 
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1 Introduction 

This paper presents a hierarchical Bayesian approach to the estimation of components' 



lO ' reliability using a Weibull model for each component in series and parallel systems. A 



series system is a frame of components that works if and only if all its components are 
functional, that is, whenever one fails the system fails. As a dual frame, the parallel 



(N . 

O . system fails if and only if all components are malfunctioning 
^ I The literature dealing with the problem of estimating the reliability of series systems, 

or competing risks, is abundant. The work of Kaplan and Meier [sl is arguably the most 
celebrated work, where it was developed a nonparametric estimator using a frequentist ap- 
^ ■ proach. For a Bayesian counterpart, we draw the attention to Salinas- Torres et al. ^] and 
o3 . Polpo and Sinha jsj. Rodrigues et al. jsj performed a simulation study of three different 
methods to estimate the reliability of a series system. They compared the Kaplan-Meier 
estimator jsf, maximum likelihood estimator (MLE) and the Bayesian plug-in estimator 
(BPE) for Weibull reliability systems. Their results indicated that MLE and BPE are 
similar in quality and that both outperformed the Kaplan-Meier estimator. However, the 
construction of credible bounds was not addressed in their work. 

For parallel systems, the literature is scarce. To the best of our knowledge, Polpo and 
Pereira jsj were the first to address the nonparametric reliability estimation in parallel 
systems and their components, using the Bayesian paradigm. Later on, Bhattacharya 
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and Samaniego [jj proposed a frequentist nonparametric estimator for components' reli- 
ability under a restrictive condition that all components are independent and identically 
distributed. 

In a related work, Polpo et al. |3] presented the reliability estimation with Weibull 
models and non-informative priors, also using the Bayesian paradigm. Their proposal 
had very demanding computational needs, and the described Markov Chain Monte Carlo 
(MCMC) algorithm suffered from convergence issues in many problem instances, making 
the estimation a difficult task. The authors realized later on that it is often impossible 
to evaluate the components' reliability with that method because of such issues, and 
unfortunately one cannot predict when the algorithm will succeed. The present work 
aims at achieving a robust estimation procedure for the estimation problem. We suggest 
to consider non-informative priors in an one-level hierarchical model as means to solve 
the estimation issues faced by the algorithm of Polpo et al. 0]. Moreover, an important 
goal of this work is to provide a simple way to build credible bounds for the reliability 
functions, giving a step-by-step algorithm that performs such analysis. 



Similarly to Polpo et al. 17] and Rodrigues et al. p\, we consider Weibull statistical 
models for the system's reliability. Two independent gamma distributions are considered 
in the first hierarchical level. The gamma distributions are parametrized by their means 
and variances, instead of the more common parametrization with scale and shape. In the 
second level of the hierarchy, we choose two flat priors for the means, and fixed values 
for the variances of the gamma distributions in the first level. These hyper-parameters 
corresponding to variances in the first level can be seen as prior precision parameters. 
The means of the gamma distributions (first level of the hierarchy) can be viewed as the 
prior expectations of the parameters of the Weibull model. In this work, the posterior 
modes are taken as the Bayesian estimators of the gamma distributions. 

The estimation of the reliability functions has three main steps: (i) we draw a sam- 
ple from the posterior distribution of the Weibull parameters; (ii) using the appropriate 
transformation, we build a sample from the reliability posterior distribution; and (iii) lo- 
cally, for each reliability time, we evaluate the posterior mean. The high posterior density 
(HPD) procedure was used to define the credible region for the reliability function. We 
emphasize that we are not using the plug-in estimator, but the posterior mean of the 
reliability function, which seems more suitable under the Bayesian paradigm. 

This paper is organized as follows. Section [2] describes all functions that are involved 
in the estimation procedure. Section [3] provides the estimation procedure itself. Section 
m presents a simulation study that highlights the quality of the model and the proposed 
estimators, and final remarks and additional comments are given in Section [51 We note 
that an extended abstract of this work has appeared in the Brazilian Conference on 
Bayesian Statistics 



We use the same notation as in Polpo et al. [7j. Consider a system of k components and 
let Xj, j G {1, . . . , A;}, be the sequence of failure times of all components. We assume that 
this sequence is composed of independent random variables with a (possibly) different 
Weibull distribution for each component. Recall that we only observe a random vector 
of two variables, namely, {T,S) with T = min(Xi, . . . , X^) for the series system and 
T = max(Xi, . . . , Xk) for the parallel system, with S = j if T = Xj, for j = 1, . . . ,k. 
The 6 quantity can be viewed as an indicator function of the component that caused the 



2 The model 
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system failure. 

Consider a sample of n independent and identically distributed systems (either all 
series or all parallel systems). The observations are represented by (T, 5) = {(Tj,5j) : i = 
1, . . . ,n}. The reliability function of the jth component is given by Rj{t) = P{Xj > t), 
j — 1, . . . ,k. Therefore, the rehabihty function is R{t) = Y['j=i ^ji^) for the series system 

and R{t) = 1 — 11^=1(1 ~ ^ji't)) for the parallel system. 

We define random variables Xj for the components' reliability with WeibuU distribu- 
tions parametrized by 9j — {Pj,r]j), that is, 



P{Xj > x\ej) = R{x\ej) = exp 



for X > 0, f3j > (shape) and r]j > (scale). Then, the likelihood function for the series 
system is given by 

k n 



and for the parallel system. 



j=i i=i 



k n 



(3) 



j=i i=i 

where / is the density function of a random variable with WeibuU distribution, — 
{6i, . . . , 6k), and 1^ is the indicator function of the set A. 

The prior distributions were considered independent with /3j ~ gamma(m/3., v^.), r]j ~ 
gamma (m,,^., f,,^.), 7r(m^J oc 7i{mrij) oc 1, and vp^ and as known constants, j = 1, . . . , k. 
Then 

n{'d) (X n{6\nif3,Vp,rnr,,Vr,)7T{rni3)7i{v/3)n{rnr,)n{Vr,) 

k 

oc Yl TT{9j\m^. , vp^ , . , 'i;^j7r(m;3 . )7r(^;^ . )7r(m^ . )7r(v^ . ) 



oc 



k . 



n 



expf-m^./^j/v^jTy- exp{-m^.77j/'f;^.} 



where "d = {6, m^, v^, m,,, i;^), = (m^^, . . . , mp^,) and = (m^^, . . . , m,,^) are the 
prior mean parameters, — {v^^, . . . , v^^,) and = (1)^^, . . . ,Vr,f.) are the prior variance 
(precision) parameters, and mp.^mr^.^vii.^Vr,^ > 0, j = 1, . . . , k. 

In this case, we have that the posterior distributions of series and parallel systems are, 
respectively. 



7r(t?|t, 6) oc7r(i?)]^ 



1=1 



Vj 



exp<- 



exp < — 



and 



7r{'d\t,S) oc 7r(t?) JJ 



i=l 



Vj 



exp < - 



1 — exp < — 



where t = (ti, . . . , t„) are the observed failure time of the system, d = (5i, . . . , 5„) are the 
indicators of which component failed, and the other quantities are as defined before. 
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3 Estimation 

For the estimation, we use the EM algorithm jil to obtain the posterior mode as estimates 
of and m^, and the MCMC procedure to generate a sample from the posterior 
distribution of f3 and rj. The algorithm steps for the estimation are briefly given as 
follows: 

1. Choose the prior precision values and v^. Note that one can set to the same 
value all the precision values, that is, f = i7/3 = v^. We suggest the use of f = 4. 

2. Choose the initial guess for the parameters to be estimated: m/3, rrir,, /3 and 77. 

3. Using the initial guess, consider m/3 and as fixed values. Employing the MCMC 
tool, generate a sample (of size Up) from the posterior distribution of (3 and 77. We 
suggest the use of Up = 1000. It may also be necessary to use a "burn-in" and a 
"jump" to ensure convergence of the MCMC. 

4. (Expectation step of the EM) Using the posterior sample of (3 and r] obtained in 
Step El evaluate the mean of the likelihood function, obtaining a "mean" function 
of and nirj. 

5. (Maximization step of the EM) Find the values of and that maximize the 
function obtained in Step SI 

6. Update the initial guess of m^j and rrirf with the values obtained in Step[5l and the 
values of (3 and rj with their posterior mean obtained in Step [31 

7. Repeat Steps |3l IH [51 and El until convergence of m/j and m,, is reached. We suggest 
the use of a tolerance value of three decimal places between the previous values of 

and and the current ones, in order to decide whether to stop iterating. 

8. Once convergence of mp and is reached, use their values to generate a sample 
from the posterior distribution of f3 and rj by applying the MCMC tool. 

Using this algorithm, we obtain the posterior mode (the Bayesian estimate) of 
and m^, and a sample (of size Hp) from the joint posterior distribution of f3 and rj. If we 
estimate the reliability function of any component (let us arbitrarily choose component 1), 
then we can notice that the estimations of the other components' reliability functions are 
very similar and could be omitted here. Consider that the sample from the posterior of the 
parameters of component 1 can be expressed as (/3ii, /3i2, • • • , f^iup) and {rju, 7712, . . . , rjin^). 
To obtain the reliability estimates and credible regions, consider the functions Yi{t) = 
F(t I f3ie,riii), £ = l,...,np, where F(t \ Pu^Vu) are Weibull distribution functions 
conditioned on (5u and rju. Consequently, the posterior mean estimate of the component's 
reliability can be expressed as 



Hence, for each fixed t, Y{t) = (Yi(t), . . . , Ynp{t)) is a sample from the posterior of the 
component I's distribution function and, to obtain the credible region, we can either use 
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the quantiles of Y{t) or evaluate the high posterior density credible interval. To estimate 
the mean reliability time, we have 

E[T I data] = — ^ E[T | (3u, r]u], i= l,...,np, 

1=1 

where E[T | Pie,rji£\ = rju^{l + 1/Pu) is the mean of a random variable with Weibull 
distribution, and r(-) is the gamma function. Note that similar procedures can be used 
to evaluate other quantities of interest; as an example, the posterior median could be 
evaluated. 

4 Examples 

Example 1 

Consider three random variables Xi, X2, X3 such that Xi has Weibull distribution with 
mean 2 and variance 4, X2 has gamma distribution with mean 2 and variance 0.667, and 
X3 has log-normal distribution with mean 2.014 and variance 6.968. We have generated 
a sample (with size n = 100) of series systems with these three components and another 
sample (again with size n = 100) of parallel systems with the exactly same three compo- 
nents. The components were chosen in order to have similar means but different variances 
and, consequently, different distributions. We have used the same theoretical components 
in both simulations (series and parallel systems) to verify, in each situation, the differ- 
ences that are due to the distinct system models with the available data. The simulated 
data have the following characteristics: (i) for the series systems, we have obtained 64%, 
80%, and 56% of censured data for components 1, 2, and 3, respectively; and (ii) for the 
parallel systems, we have observed 61%, 68% and 71%, respectively for the same three 
components. In this case, the main interest is in the estimation of the components' re- 
liability functions. Note that, with our simulated example, we have a huge amount of 
censored data, making it a challenging example. 

As already said, the estimation procedures are performed using MCMC We have 
discarded the first 10, 000 samples (as burn-in) from the posterior to achieve the sta- 
tionary measure and then have generated a sample from the posterior. To perform the 
estimation of the reliability functions and the credible region, we have used a sample of 
size 1,000 from the posterior, which was obtained by discarding 10 samples (the jump 
between each final sample point). We have used f = 4 in the prior specification for 
all parameters and systems. For the experiment with series systems, we have obtained 
= (1.31,4.12,1.44) and m;, = (2.19,2.03,1.80). For the parallel systems, the esti- 
mates are = (1.28,2.57,0.92) and = (2.67,2.28,2.11). To evaluate the quality of 
these estimates, we have compared the "true" reliability of each component with the esti- 
mated reliability function. Table [T] presents the posterior mean and the posterior standard 
deviation of each parameter involved in the model (for both series and parallel systems). 
Note that the standard deviations are relatively small, and the estimation of the mean 
reliability time is very close to the original values, indicating a good performance of our 
method. 
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Table 1: Posterior mean (sd) of some quantities involved in the estimation of the simulated 
examples. 



Series system estimates 
/3 T] E{T\(3,r,) 



Component 1 
Component 2 
Component 3 



1.26 (0.17) 
3.98 (0.53) 
1.40 (0.17) 



2.27 (0.41) 
2.06 (0.12) 
1.83 (0.22) 



2.13 (0.45) 
1.87 (0.11) 
1.68 (0.22) 



Parallel system estimates 
13 r, E(T 1/3,77) 



Component 1 
Component 2 
Component 3 



1.23 (0.15) 
2.47 (0.32) 
0.87 (0.11) 



2.60 (0.27) 
2.25 (0.14) 
1.98 (0.33) 



2.45 (0.22) 
2.00 (0.13) 
2.17 (0.29) 



It can be seen from the 95% credible bounds that the "true" reliability of each com- 
ponent was well estimated. We note however that the "true" reliability functions of the 
components are, for short (time) intervals, outside the 95% credible bounds. Considering 
that these are reasonably challenging examples, this situation is likely to happen in any 
estimation procedure (see Figures [T] and [2]) . 






Figure 1: Reliability of the components in the experiment with series systems: (a) Com- 
ponent 1; (b) Component 2; (c) Component 3. 



Component 1 Component 2 Component 3 




0369 12 02468 0369 12 



Figure 2: Reliability of the components in the experiment with parallel systems: 
Component 1; (b) Component 2; (c) Component 3. 



Example 2 

This example has a simulation study to show the quality of the proposed hierarchical 
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Weibull model over many different conditions. We have considered 108 different scenar- 
ios that were built using three different sample sizes (n = 30, 100, 1000), three different 
proportions of censored data (0%, 20% and 40%), two different means (2 and 7) of reliabil- 
ity time for the generating distribution, three different generating distributions (Weibull, 
gamma, log-normal), and two types of censoring (right and left). In all these cases, we 
have fixed the variance of the generating distribution to 5. We have used a non-random 
censored approach to guarantee the desired proportions of censure, that is, for each sce- 
nario, we have fixed a time for which all values that are larger than this time are assumed 
censored for the right-censored data (series systems), and all values smaller than this 
fixed time are assumed censored for the left-censored data (parallel systems). In order to 
improve the analysis, 100 copies of each scenario were considered. 

To summarize and to compare the results of the simulation, we have evaluated the 
bias and the mean squared error (MSE) of the estimated mean reliability time, for each 
scenario. The results are presented in Tables [2J171 One can see from the results that 
all biases and MSEs are close to zero, indicating that, in all scenarios, the model has 
estimated well the true mean reliability time. Even in the most challenging scenarios, 
which are those with small sample size [n = 30) and large proportion of censoring (40%), 
the estimated mean reliability time has had a good performance. When the generating 
distribution is Weibull, the model should obviously estimate well, yet the results show that 
the performances were good even for the other two generating (non- Weibull) distributions. 

Table 2: Bias and mean squared error (MSE) of the E(T) estimate for data generated 
from the Weibull distribution with right-censored data. 







n = 


30 


n = 


100 


n = 


1000 


Gens. 


True E(T) 


Bias 


MSE 


Bias 


MSE 


Bias 


MSE 


0% 


2 


-0.1300 


0.2217 


-0.0426 


0.0572 


-0.0057 


0.0043 


0% 


7 


0.0293 


0.1523 


0.0026 


0.0474 


-0.0095 


0.0043 


20% 


2 


-0.2813 


0.4893 


-0.1242 


0.1275 


-0.0200 


0.0068 


20% 


7 


-0.0522 


0.2307 


-0.0233 


0.0565 


-0.0148 


0.0047 


40% 


2 


-0.4293 


0.7364 


-0.2394 


0.2925 


-0.0239 


0.0108 


40% 


7 


-0.1773 


0.3590 


-0.0661 


0.0874 


-0.0162 


0.0055 



Table 3: Bias and mean squared error (MSE) of the E(T) estimate for data generated 
from the gamma distribution with right-censored data. 







n = 


30 


n = 


100 


n = 


1000 


Gens. 


True E(T) 


Bias 


MSE 


Bias 


MSE 


Bias 


MSE 


0% 


2 


-0.1407 


0.2286 


-0.0191 


0.0593 


-0.0003 


0.0055 


0% 


7 


0.0006 


0.1845 


0.0081 


0.0623 


0.0098 


0.0048 


20% 


2 


-0.3829 


0.6346 


-0.1132 


0.1234 


-0.0560 


0.0111 


20% 


7 


0.0569 


0.1721 


0.1121 


0.0669 


0.1049 


0.0153 


40% 


2 


-0.5068 


0.9424 


-0.3084 


0.3807 


-0.1174 


0.0299 


40% 


7 


0.1389 


0.2505 


0.2372 


0.1107 


0.2336 


0.0582 
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Table 4: Bias and mean squared error (MSE) of the E(T) estimate for data generated 
from the log-uoiiiial distribution witli riglit-ceusorcxl data. 







n — 


'in 
oU 


n — 


1 nn 

iUU 


n — 


1 nnn 

iUUU 


Cens. 


True E(r) 


Bias 


MSE 


Bias 


MSE 


Bias 


MSE 


0% 


2 


-0.1160 


0.1757 


-0.0024 


0.0396 


-0.0089 


0.0050 


0% 


7 


-0.0348 


0.1551 


0.0462 


0.0429 


0.0220 


0.0052 


20% 


2 


0.0720 


0.1410 


0.2475 


0.0779 


0.2523 


0.0655 


20% 


7 


0.0786 


0.1492 


0.1848 


0.0637 


0.1620 


0.0296 


40% 


2 


0.2373 


0.1795 


0.4205 


0.1981 


0.4371 


0.1928 


40% 


7 


0.2190 


0.2094 


0.3438 


0.1532 


0.3337 


0.1147 



Table 5: Bias and mean sqiiarcd error (MSE) of the E(T) estimate for data generated 
from the Wcibull distribution witli Icft-ccnsorcd data. 







n = 


30 


n = 


100 


n = 


1000 


Cens. 


True E(r) 


Bias 


MSE 


Bias 


MSE 


Bias 


MSE 


0% 


2 


-0.1300 


0.2217 


-0.0426 


0.0572 


-0.0057 


0.0043 


0% 


7 


0.0293 


0.1523 


0.0026 


0.0474 


-0.0095 


0.0043 


20% 


2 


0.1030 


0.1178 


0.0034 


0.0475 


-0.0060 


0.0043 


20% 


7 


0.0387 


0.1562 


-0.0001 


0.0490 


-0.0094 


0.0043 


40% 


2 


-0.1813 


0.2660 


-0.0471 


0.0579 


-0.0065 


0.0043 


40% 


7 


0.0301 


0.1813 


0.0040 


0.0597 


-0.0099 


0.0048 



Table 6: Bias and mean squared error (MSE) of the E(T) estimate for data generated 
from the gamma distribution with left-censored data. 







n 


30 


n = 


100 


n = 


1000 


Cens. 


True E(T) 


! ') ■ . 


MSE 


Bias 


MSE 


Bias 


MSE 


0% 


2 


-0.1407 


0.2286 


-0.0192 


0.0594 


-0.0003 


0.0055 


0% 


7 


0.0006 


0.1845 


0.0081 


0.0623 


0.0098 


0.0048 


20% 


2 


0.1029 


0.1038 


0.0448 


0.0490 


0.0008 


0.0054 


20% 


7 


0.1108 


0.2143 


0.0950 


0.0739 


0.0932 


0.0139 


40% 


2 


-0.1934 


0.2463 


-0.0245 


0.0602 


-0.0029 


0.0055 


40% 


7 


0.2185 


0.2842 


0.2105 


0.1242 


0.1850 


0.0402 



Table 7: Bias and mean squared error (MSE) of the E(T) estimate for data generated 
from the log-normal distribution with left-censored data. 







n 


30 


n 


100 


n = 


1000 


Cens. 


True E(r) 


1 


MSE 


1 


MSE 


Bias 


MSE 


0% 


2 


-0.1160 


0.1757 


-0.0024 


0.0396 


-0.0089 


0.0050 


0% 


7 


-0.0348 


0.1551 


0.0462 


0.0429 


0.0220 


0.0052 


20% 


2 


-0.0092 


0.1253 


0.0254 


0.0405 


0.0211 


0.0052 


20% 


7 


0.1003 


0.1701 


0.1757 


0.0784 


0.1375 


0.0240 


40% 


2 


-0.1165 


0.1937 


0.0575 


0.0441 


0.0624 


0.0086 


40% 


7 


0.2359 


0.2727 


0.3138 


0.1548 


0.2763 


0.0825 
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5 Final Remarks 

We have introduced a Bayesian reliability statistical analysis using hierarchical models 
for the problem of estimating the reliability functions and credible bounds of series and 
parallel systems. The MCMC has shown good performance in terms of convergence, 
making the inference process simple and efficient. It shall be noted that this performance 
is not dependent on our choice of a "non-informative" scheme to define the prior hyper- 
parameters. This is important because other researchers may want to fairly compare 
our method with other frequentist estimators. However, informative priors may very 
well produce additional improvements in the estimates. The Example [1] has shown good 
robustness in the sense that the model has performed well for all components in both series 
and parallel systems. Another important aspect is that we can obtain credible bounds 
for the reliability function, task that is usually hard if one uses a plug-in estimator for 
the reliability function. The Example [2] provides an extensive simulation study with more 
than one hundred different scenarios. Overall, the model has performed very well for 
estimating the mean reliability time. Some open questions that should be addressed in 
future works are the development of hypothesis tests for the components, for instance, 
one can have interest in testing the hypothesis of equal means of all components (or a 
subset of components), and the extension of these ideas to more general coherent systems. 
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